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Abstract 

Quantum  spin  models  with  spatially  dependent  interactions,  known  as  compass  models,  play  an 
important  role  in  the  study  of  frustrated  quantum  magnetism.  One  example  is  the  Kitaev  model  on  the 
honeycomb  lattice  with  spin-liquid  (SL)  ground  states  and  anyonic  excitations.  Another  example  is 
the  geometrically  frustrated  quantum  120°  model  on  the  same  lattice  whose  ground  state  has  not  been 
unambiguously  established.  To  generalize  the  Kitaev  model  beyond  the  exactly  solvable  limit  and 
connect  it  with  other  compass  models,  we  propose  a  new  model,  dubbed  The  tripod  model’,  which 
contains  a  continuum  of  compass-type  models.  It  smoothly  interpolates  the  Ising  model,  the  Kitaev 
model,  and  the  quantum  120°  model  by  tuning  a  single  parameter  6\  the  angle  between  the  three  legs 
of  a  tripod  in  the  spin  space.  Hence  it  not  only  unifies  three  paradigmatic  spin  models,  but  also  enables 
the  study  of  their  quantum  phase  transitions.  We  obtain  the  phase  diagram  of  the  tripod  model 
numerically  by  tensor  networks  in  the  thermodynamic  limit.  We  show  that  the  ground  state  of  the 
quantum  120°  model  has  long-range  dimer  order.  Moreover,  we  find  an  extended  spin-disordered 
(SL)  phase  between  the  dimer  phase  and  an  antiferromagnetic  phase.  The  unification  and  solution  of  a 
continuum  of  frustrated  spin  models  as  outline  here  maybe  useful  to  exploring  new  domains  of  other 
quantum  spin  or  orbital  models. 


1.  Introduction 

Model  Hamiltonians  describing  interacting  spins  localized  on  lattice  sites  are  at  the  central  stage  in  the  field  of 
quantum  magnetism.  A  class  of  spin  models,  collectively  known  as  the  compass  models  [1],  stand  out  owing  to  a 
unique  feature  they  share  in  common:  the  spin  exchange  interactions  differ  for  different  lattice  bond 
orientations.  This  is  in  contrast  to  the  familiar  Heisenberg  model  or  the  Ising  model,  where  the  exchange  has  the 
same  form  for  all  bonds  connecting  the  nearest  neighboring  sites.  The  compass  models  arise  naturally  as  low 
energy  effective  Hamiltonians  in  Mott  insulators  with  orbital  degrees  of  freedom  [2-7]  as  well  as  interacting 
systems  with  spin-orbit  coupling.  These  highly  nontrivial  models  are  also  very  appealing  from  a  pure  theoretical 
point  of  view  because  they  offer  a  natural  arena  to  study  frustrated  quantum  magnetism  [8, 9] .  Exactly  solvable 
compass  models,  the  Kitaev  model  in  particular,  have  played  a  pivotal  role  in  stimulating  the  field  of  topological 
quantum  computing  [10, 1 1].  The  rich  physics  contained  in  compass  models  has  been  reviewed  recently  in  [1]. 

Our  work  is  directly  motivated  by  two  well  known  compass  models  defined  on  the  honeycomb  lattice.  The 
first  example  is  the  Kitaev  model  [11],  where  the  exchange  interactions  between  two  neighboring  sites  are  given 
by  a  f  a* ,  o\  cryp  and  of  crz  respectively.  As  shown  by  Kitaev,  this  model  is  exactly  solvable  and  has  anyonic 
excitations  obeying  fractional  statistics  [11].  The  spatially  dependent  exchange  interactions  suppress  long-range 
spin  order  and  support  a  quantum  spin  liquid  (SL)  ground  state,  one  of  the  most  sought  after  exotic  many-body 
states  in  condensed  matter  physics  [12].  The  Kitaev  model,  despite  its  theoretical  appeal,  is  neither  readily 
realized  in  materials  nor  easily  simulated  with  synthetical  quantum  matter  such  as  cold  atoms  on  optical  lattices. 
Recently,  the  hybrid  Kitaev-Heisenberg  model,  a  linear  superposition  of  a  Kitaev  term  and  a  Heisenberg  term, 
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was  proposed  for  iridium  oxides  and  solved  numerically  [13].  Besides  the  SL  phase,  the  hybrid  model  contains 
other  interesting  phases  such  as  the  stripe  and  the  zigzag  phase  due  to  the  competition  between  the  two  terms 
[13].  The  phase  diagram  becomes  even  richer  when  off-diagonal  spin  exchange  interactions  are  added  [14-16]. 

The  second  example  of  compass  models  is  the  quantum  120°  model  [6,  ?].  It  is  very  analogous  to  the  Kitaev 
model  but  the  spin  operators  <jx ,  <jy ,  and  az  for  the  three  bond  directions  are  replaced  by  three  (pseudo)spin  1/2 
operators  T1,  T2,  and  T3  respectively.  They  form  an  angle  of  120°  with  each  other  on  the  xz  plane  in  spin  space, 
hence  the  name  The  120°  model’.  It  was  introduced  to  described  the  low  energy  physics  of  transition  metal 
oxides  [17]  with  doubly  degenerate  orbitals,  e.g.  orbital- only  models  of  eg  orbitals  on  cubic  lattice  [3].  Later,  two 
of  us,  and  Wu  independently,  found  that  the  120°  model  can  be  naturally  realized  in  strongly  interacting  spinless 
p-orbital  fermions  on  the  honeycomb  optical  lattice  [6, 7].  Although  it  is  geometrically  frustrated,  spin  wave 
analysis  indicates  that  long-range  order  is  stabilized  by  quantum  fluctuations  through  the  order  by  disorder 
mechanism  [6,  7].  While  the  semiclassical  spin- wave  analysis  is  suggestive,  the  ground  state  of  the  120°  model  on 
honeycomb  lattice  remains  to  be  identified  unambiguously. 

Given  the  apparent  similarities  between  the  Kitaev  model  and  the  120°  model,  it  is  natural  to  seek  the 
conceptual  and  quantitative  link  between  them.  Indeed,  these  two  models  can  be  viewed  as  two  instances  of  a 
more  general  class  of  compass  models  [18].  In  this  paper,  we  provide  a  concrete  construction  and  propose  a 
‘super  model’  which  contains  three  paradigm  models,  the  Ising,  the  Kitaev  and  the  120°  model,  as  special  limits. 
It  only  has  a  single  tuning  parameter  O'  and  a  simple,  intuitive  picture  for  the  three  (pseudo)spin  operators:  they 
form  a  tripod  in  spin  space  as  shown  in  figure  1 .  Analogous  to  tuning  the  tripod  to  raise  or  low  a  mounted  camera 
in  photography,  dialing  the  angle  O'  between  the  three  legs  takes  the  Ising  model  (tripod  fully  closed)  smoothly  to 
the  Kitaev  model  (tripod  open  with  three  legs  orthogonal  to  each  other)  and  then  to  the  120°  model  (tripod  fully 
open  with  three  legs  in  the  same  plane).  Immediately,  one  conjectures  that  the  phase  diagram  of  this  continuum 
compass  model  is  highly  nontrivial  containing  drastically  different  long-range  ordered  states  as  well  as  SLs. 

We  obtain  the  phase  diagram  of  this  ‘super  model’  using  tensor  network  states  which  have  gained 
considerable  success  recently  in  the  study  of  frustrated  magnetism  [19-22].  The  results  are  summarized  in 
figures  1  and  2.  The  order  parameters  are  calculated  using  the  tensor  renormalization  group  (TRG)  method 
formulated  in  thermodynamic  limit  [23, 24].  We  show  that  the  ground  state  of  the  quantum  120°  model  is  a 
long-range  ordered  dimer  phase,  and  a  SL  phase  exists  in  an  extended  region  in  our  phase  diagram  .  The 
numerical  results  of  TRG  are  further  confirmed  and  crosschecked  with  projected  entangled  pair  states  (PEPS) 
calculations  [25, 26]  for  finite  systems,  exact  diagonalization  (ED),  and  spin  wave  analysis.  We  discuss  the 
qualitative  features  of  the  quantum  phase  transitions  between  the  SL  phase  and  the  dimer  phase  by  introducing  a 
topological  charge  (spin  vortex)  for  the  dimer  configuration.  We  further  show  that  the  proposed  tripod  model 
can  in  principle  be  simulated  with  Hubbard  model  in  the  Mott  insulating  regime,  e.g.,  using  cold  atoms  on 
optical  lattice  with  artificial  gauge  fields. 

2.  The  tripod  model 

We  generalize  the  Kitaev  and  the  120°  model  to  the  following  continuum  compass  model  defined  on  the  two- 
dimensional  honeycomb  lattice 

H(d)  =  jY:V(0)V+eye)  d) 

r,7 


where  /  >  0,  and  for  each  lattice  site  r  the  spin  1/2  operators  are  defined  as 


T7(6>) 


—  ( rz  cos  6  +  rx  sin  6  )  cos  0  +  — ry  sin  0 
2  7  7  2 


(2) 


with  the  Pauli  matrices  rx,  ry ,  rz.  Each  site  r  is  coupled  to  its  neighbors  r  +  e7,  where  e7,  7  =  1,  2,  3,  denotes 
the  three  bond  vectors  of  the  honeycomb  lattice.  Geometrically,  the  three  T7  form  a  tripod  in  the  spin  space  as 
shown  in  figure  1 :  they  are  tilted  from  the  xz  plane  by  angle  0  and,  when  projected  onto  the  xz  plane,  are 
orientated  along  the  corresponding  bond  direction  e7,  i.e.  at  azimuthal  angle  07  =»  0,  27 r /3,  47 r /3  respectively. 
While  T7  is  most  naturally  defined  through  the  tilting  angle  0,  it  is  much  more  convenient  to  introduce  another 
angle,  O',  to  discuss  the  various  limits  of  H  ( 0 ).  O'  is  the  angle  between  T1  and  T2,  i.e.,  the  two  adjacent  legs  of  the 
tripod.  And  it  is  related  to  0  by  trigonometry 

5  In  this  paper,  we  use  the  term  ‘spin  liquid’  to  denote  a  phase  that  shows  no  long-range  spin  order  according  to  our  tensor  network 
algorithms.  At  the  special  point,  9'  =  90°  the  tripod  model  reduces  to  the  Kitaev  model  and  its  ground  state  is  well  established  to  be  a  gapless 
spin  liquid.  Our  numerical  results  offer  evidence  that  the  same  spin-liquid  state  will  survive  away  from  the  Kitaev  point.  The  nature  of  the 
excitations  (e.g.  their  fractional  statistics)  remains  to  be  checked  in  order  to  firmly  establish  that  it  is  a  quantum  spin  liquid. 
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Figure  1.  The  tripod  model  on  the  honeycomb  lattice,  equation  (1).  (a)  The  nearest  neighbor  spin  exchange  along  bond  direction 
e7,  7=  1,  2,  3,  is  defined  through  the  spin  1/2  operator  T7,  represented  by  an  arrow  in  spin  space  spanned  by  rx,  ry ,  rTThethree 
T7  can  be  thought  as  the  three  legs  of  a  tripod,  being  tilted  out  of  the  xz  plane  by  angle  9 ,  and  forming  an  angle  9'  with  each  other. 
When  projected  onto  the  xz  plane,  T7  is  along  the  e7  direction,  (b)  The  schematic  phase  diagram  of  the  tripod  model.  As  the  tripod  is 
opened  by  increasing  9',  the  model  starts  as  the  Ising  model  at  9'  =  0,  becomes  the  Kitaev  model  at  9'  =  90°,  and  then  the  quantum 
120°  model  at  9'  =  120°.  Three  phases  are  identified:  a  Neel  ordered  antiferromagnet,  a  spin  liquid  (SL),  and  a  long-range  ordered 
dimer  phase. 


Figure  2.  Identifying  the  phases  of  the  tripod  model.  The  plot  shows  the  order  parameters  Ch  (filled  squares)  and  02  (empty  circles)  as 
a  function  of  9'  calculated  from  the  infinite  tensor  network  algorithms  with  bond  dimension  D  =  8  and  6-sites  unit  cell.  The  insets 
illustrate  the  schematic  spin  configurations  in  the  Neel  ordered  phase  (left)  and  the  dimer  phase  (right)  for  a  hexagon.  For 
9'  G  [87°,  94°],  the  spin  averages  are  zero,  and  the  ground  state  is  identified  as  a  spin  liquid. 
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cos  O'  =  1 - cos2 

2 


(3) 


We  will  take  O'  as  the  only  tuning  parameter  in  the  tripod  model. 

Three  special  limits  of  this  model  can  now  be  identified.  First  of  all,  when  O'  =  0,  T7  all  collapse  to  ry.  The 
tripod  is  closed,  and  H  (0)  is  nothing  but  the  Ising  model 


4  r,7 


(4) 


Note  that  we  choose  ry  as  the  vertical  axis  in  spin  space  instead  of  the  usual  convention  of  rz  so  that  H  ( 0 ) 
reduces  exactly  to  the  120°  model  defined  in  our  earlier  work  [6].  Secondly,  when  O'  =  90°,  H  ( 0 )  reduces  to  the 
Kitaev  model  since  the  operators  T7  are  now  perpendicular  to  each  other  in  the  spin  space.  We  can  simply 
identify  them  as  ax,  ay  and  az  (apart  from  a  factor  1  /2)  in  a  rotated  coordinate  system  as  illustrated  in 
figure  1(b).  Thirdly,  for  O'  =  120°,  H  (0)  becomes  the  quantum  120°  with  T1,  T2,  T3  all  confined  within  thexz 
plane.  It  can  be  visualized  as  a  fully  open  tripod. 

As  well  known,  the  Ising  model  has  antiferromangetic  (AF)  order  with  the  order  parameter 


Oi  =  (It'D/ 2. 


(5) 


On  the  other  hand  side,  the  quantum  120°  model  is  conjectured  to  be  long-range  ordered  despite  the  geometric 
frustration.  We  introduce  the  following  ‘order  parameter’  to  measure  the  in-plane  magnetization 

02  =  V(r*)2  +  (tz)2  /2.  (6) 


By  solving  H  (0)  using  tensor  network  algorithms,  we  compute  the  average  spin  (r^,y,z)  in  the  ground  state. 
The  main  results  are  summarized  in  the  schematic  phase  diagram  in  figure  1  (b).  Figure  2  shows  the  variation  of 
the  two  order  parameters  introduced  above  as  O'  is  changed.  The  region  at  small  O'  corresponds  to  the  familiar 
Neel  order  which  is  characteristic  of  the  classical  Ising  model  and  illustrated  in  the  left  inset  of  figure  2.  Despite 
the  increased  quantum  fluctuations  as  O'  is  increased,  the  Neel  ordered  phase  persists  up  to  O'  ~  87°.  At  the 
opposite  end  of  large  0' ,  we  find  that  the  long-range  spin  order  consists  of  a  set  of  £dimers,’  i.e.  opposite  spins  on 
neighboring  sites,  arranged  into  a  periodic  pattern  of  triangular  lattice  (figure  5(a)).  The  triangular  lattice  and  its 
enlarged  unit  cell  becomes  transparent  if  we  introduce  a  topological  charge  (red  dot  in  figure  5(a))  for  each 
hexagon  with  spins  all  pointing  outwards.  If  we  focus  on  one  individual  hexagon,  e.g.  the  one  shown  in  the  right 
inset  of  figure  2,  the  orientations  of  the  dimers  happen  to  be  also  60°  (or  equivalently  120°)  apart.  We  will  refer  to 
this  phase  simply  as  the  £dimer  phase.’  In  particular,  it  is  the  ground  state  of  the  quantum  120°  model  on 
honeycomb  lattice.  This  point  will  be  further  discussed  in  section  5. 

Sandwiched  between  the  Neel  ordered  phase  and  the  dimer  phase,  a  quantum  SL  phase  is  stabilized  for 
O'  G  [87°,  94°].  The  conclusion  is  mainly  based  on  the  observation  from  figure  2  that  the  order  parameters  Oy2 
in  this  region  are  nearly  zero  compared  to  those  in  other  two  phases.  This  conclusion  is  also  consistent  with  the 
exactly  solution  of  Kitaev  model  for  O'  =  90°.  The  order  parameters  as  functions  of  O'  in  figure  2  also  suggest 
that  the  two  quantum  phase  transitions  in  the  tripod  model  maybe  qualitatively  different.  The  gradual  drop  of 
02  at  0'  >  90°  indicates  a  continuous  phase  transition  between  the  dimer  phase  and  the  SL  phase.  In  contrast, 
the  drop  of  the  order  parameter  Ox  at  0'  <  90°  is  rather  sharp,  pointing  to  a  likely  first-order  phase  transition. 
The  details  of  the  calculations  leading  to  these  results  will  be  discussed  below  in  section  3. 


3.  Tensor  network  algorithms 

Recent  developments  of  entanglement-based  tensor  network  algorithms  provide  a  novel,  accurate  approach  to 
strongly  correlated  electron  systems  [i  —29].  Particularly,  they  have  been  successfully  applied  to  frustrated 
quantum  magnets  [19-22]  and  the  t  —  J  model  [30, 3 1]  to  yield  insights  previously  unattainable  from 
conventional  methods.  To  find  the  phase  diagram  of  the  proposed  tripod  model,  we  employ  two 
complementary  tensor  networks  algorithms,  one  for  finite-size  systems  and  the  other  for  infinite  systems  in  the 
thermodynamic  limit,  to  find  the  ground  state  and  the  order  parameters.  In  both  algorithms,  the  ground  state 
wave  function  is  constructed  as  a  network  of  local  tensors  defined  on  lattice  sites.  Each  tensor  has  one  physical 
index  representing  the  spin  degree  of  freedom  and  three  virtual  indices,  each  with  bond  dimension  D,  describing 
the  quantum  entanglement  with  its  three  neighboring  sites. 

We  first  apply  the  finite  PEPS  algorithm  [25-27]  to  solve  H  (6)  for  a  six-site  system  with  periodic  boundary 
conditions.  The  ground  state  energies  obtained  coincide  with  those  from  ED.  This  suggests  that  PEPS  is 
intrinsically  superior  compared  to  mean  field  theories  when  applied  to  frustrated  spin  Hamiltonians  such  as 
H  ( 0 ).  The  order  parameters  decay  to  zero  as  the  ground  state  is  approached  for  a  finite  system.  Nonetheless, 
their  decay  behaviors  are  quite  disparate  for  0'  values  in  the  Ising,  Kitaev,  and  120°  regions,  suggesting  three 
different  phases.  The  details  of  the  calculation  are  presented  in  appendix  A. 
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Figure  3.  The  HOTRG  procedure  used  to  calculate  the  expectation  value  (O)  after  the  convergence  of  the  ground  state,  (a)  The  unit 
cell  (the  center  hexagon)  consists  of  six  sites.  Associated  with  each  site  is  a  local  tensor  Tj  at  the  zth  RG  step,  (b)  The  coarse  graining 
procedure  to  construct  the  new  tensor  T/+ 1  from  the  old  tensors  T\  and  its  neighbors  Tj,  Tj,  Tj.  The  other  five  tensors  are  updated  in 
the  similar  way.  (c)  The  new  tensors  Tj+ 1  again  form  a  honeycomb  lattice. 


To  study  the  tripod  model  in  the  thermodynamic  limit,  we  first  find  the  converged  ground  state  using 
imaginary  time  evolution  and  following  the  simple  update  scheme  as  described  in  [32]  which  generalizes  the 
time-evolving  block  decimation  (TEBD)  [33]  technique  to  two-dimensions.  For  a  n-site  unit  cell,  e.g.  a  six-site 
unit  cell  shown  in  figure  3(a),  we  need  3n/2  different  bond  vectors  that  represent,  roughly  speaking,  a  mean- 
field  approximation  of  the  environment.  Using  these  bond  vectors,  the  simple  update  starts  with  n  random 
tensors  and  iterates  until  convergence  is  achieved.  At  the  end  of  the  calculation,  the  ground  state  |4/)  is 
characterize  by  n  tensors  Tj ,  j  E  [1 ,  n\.  We  then  evaluate  the  expectation  value  of  operator 
O,  (O)  =  (4/| 0|4/)  /  (4/|4/)  which  involves  the  (infinite)  product  of  tensors  7},  using  a  real  space  coarse  graining 
procedure  known  as  higher-order  TRG  (HOTRG)  [24]  schematically  shown  in  figure  3.  We  outline  the  main 
steps  here.  At  the  ith  step  of  HOTRG,  a  local  tensor,  say  T{,  is  regrouped  with  its  three  nearest  neighbor  tensors 
(T2,  Zi,  T6 )  to  form  a  new  tensor  Tjz+1.  More  generally,  for  odd  or  even  sites 


T’+1  =  ^T'T’T'T', 

S.l. 

i+1 


T, 


y^T/T'r', 

s.l. 


(7) 

(8) 


where  the  summation  is  over  the  shared  legs  (abbreviated  as  s.l.,  the  solid  lines  in  figure  3(b))  of  the  neighboring 
tensors,  i.e.  tensor  contractions.  The  new  tensors,  each  of  which  contains  four  old  tensors,  are  of  higher 
dimensions  and  truncated  to  have  the  same  dimension  as  T*  via 

Ti/e1  =  'Efi/e1UxUyUz,  (9) 

S.l. 
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where  the  three  projection  tensors  UXyy?z,  shown  in  figure  3(b),  are  obtained  as  follows.  Take  Ux  as  an  example. 
First,  fj  is  reshaped  into  matrix  Tj  with  the  row  corresponding  to  the  leg  along  the  x- direction.  Then,  a  matrix  hi 
is  obtained  by  singular  value  decomposition 

6 

YN,'T)  =  (10) 

;'=  i 

Finally,  Ux  is  obtained  by  truncating  hi  to  a  given  truncation  dimension  \  and  reshaping  it  back  to  the  tensor 
form.  The  new  tensors  T*+ 1  now  form  exactly  the  same  honeycomb  lattice  structure  as  the  old  tensors  T1  but 
represent  a  larger  system,  see  figure  3(c).  This  constitutes  a  single  RG  step.  By  iterating  the  RG  steps  many  times, 
the  converged  result  of  (O)  well  approximates  the  expectation  value  in  the  thermodynamic  limit. 

By  following  these  procedures,  we  have  calculated  the  ground  state  energy  and  the  ground  state  expectation 
values  of  the  order  parameters  Oh  02  for  different  unit  cell  sizes,  n  =  2,  4,  6,  8.  We  found  that  the  six- site  unit 
cell  gives  the  lowest  energy.  The  two-site  unit  cell  yields  results  in  agreement  with  the  six-site  unit  cell  within  the 
parameter  region  O'  <  90°.  The  four-site  and  eight-site  unit  cells,  however,  lead  to  excited  states  with 
significantly  higher  energy.  Thus,  we  conclude  that  the  six- site  unit  cell  is  the  most  reasonable  choices  for  all  the 
O'  values  in  the  ground  state  calculation.  In  practice,  one  can  safely  use  the  two-site  unit  cell  for  O'  <  90°  since  it 
is  significantly  cheaper.  The  phase  diagram  and  the  spin  configurations  in  the  ordered  phases  shown  in  figure  2 
are  obtained  by  using  the  two-site  unit  cell  for  O'  90°  and  the  six-site  unit  cell  for  O'  >  90°. 

4.  Spin  wave  analysis 

To  cross-check  the  TRG  results,  we  perform  the  standard  spin  wave  analysis  of  the  tripod  model.  It  is  important 
to  keep  in  mind  that  the  validity  of  the  spin  wave  theory,  which  can  be  viewed  as  expansion  in  series  of  1/5, 
becomes  questionable  in  the  limit  of  5  =  1/2.  Yet  the  analysis  offers  a  rough  picture  of  the  role  played  by 
geometric  frustration  and  how  the  Neel  order  and  dimer  order  get  destroyed  by  the  increased  quantum 
fluctuations.  As  we  will  show  below,  the  estimations  of  the  two  quantum  critical  points  from  the  spin  wave 
theory  turn  out  to  be  in  broad  agreement  with  the  phase  digram  predicted  by  the  tensor  network  algorithms. 

The  analysis  starts  by  partitioning  the  honeycomb  lattice  into  the  A  and  B  sublattice  and  introducing 
5r  =  dz  Tr  for  all  sites  on  the  A  (B)  sublattice.  Then  the  tripod  Hamiltonian  acquires  a  suggestive  form 

H(0)  =  J-  £  [S7(0)  -  S7+e7(0)]2  +  const.  (11) 

^  r£A,7 

Here  we  have  promoted  the  spin  1/2  operator  r/2  to  general  spin  operator  S  with  spin  quantum  number  5.  It 
follows  that  classical  ground  states  are  massively  degenerate  (except  for  the  Ising  limit).  Any  spin  configurations 
with  5r7  (0)  =  5rT|_e  (0),  i.e.,  the  projection  of  S  along  the  bond  direction  being  the  same  for  any  two  neighboring 
sites,  will  minimize  the  classical  energy.  This  is  a  well  known  feature  of  compass  models,  see  the  review  [1].  The 
special  case  of  the  classical  120°  model  on  honeycomb  lattice  was  previously  discussed  in  [7, 34].  We  will  confine 
our  spin  wave  analysis  to  the  simple  case  of  spatially  homogeneous  spin  configurations  Sr  =  S0  as  done  in  [34]. 
The  direction  of  S0  is  characterized  by  its  polar  angle  ip  measured  from  ry  and  its  azimuthal  angle  a  of  S0 
measured  from  rz  in  rx-rz  plane.  The  corresponding  classical  ground  state  energy  per  unit  cell  is 

=  —  —  (2  sin2  0  cos2  p  +  cos2  6  sin2  p).  (12) 

S2J  2 

It  is  interesting  to  note  that,  coincidently,  at  the  Kitaev  point,  0'  =  90°  which  corresponds  to 
0  =  cos-1  V2/3,£o  is  completely  flat  and  does  not  depend  on  p.  For  0'  <  90°,  E0  is  minimized  when  ip  =  Oor 
7 r,  corresponding  to  the  two  degenerate  states  with  spin  up  or  down  in  the  Neel  ordered  phase.  In  contrast,  for 
0'  >  90°,  E0  is  minimized  when  <p  =  7r/2,  i.e.,  S0  lies  within  the  rx-rz  plane.  Therefore,  the  mean  field  theory 
above  predicts  that  the  tripod  model  has  a  phase  transition  exactly  at  the  Kitaev  point. 

Applying  the  Holstein-Primakoff  transformation  [35]  to  H  ( 0 )  and  expanding  the  resulting  Hamiltonian  of 
bosons  to  order  1/5,  we  compute  the  quantum  fluctuation  correction  to  the  ground  state  energy  for  the  two 
long- ranged  ordered  states  respectively  and  find 

7  =  2EWk)|-A,  (13) 

^  k,A 

where  A  =  —  Eq/ 52/,  N  is  the  number  of  sites  within  the  A  sublattice,  the  k  summation  is  over  the  first 
Brillouin  zone  of  the  A  sublattice,  and  (k)  describe  the  spin  wave  dispersion  for  branch  A  =  1 ,  2  and  they  are 

given  by  the  eigenvalues  of  the  matrix 
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Figure  4.  The  energy  per  unit  cell  from  leading  order  spin  wave  theory  for  the  tripod  model.  Upper  panel:  E  {(j))  for  the  Neel  ordered 
phase.  At  9'  ~  75.0°,  the  location  and  number  of  minima  of  E  change,  indicating  the  transition  to  a  different  phase.  Lower  panel: 
E{(f)  =  7t/2,  a)  for  the  dimer  phase  with  minima  occurring  at  an  =  mr/3.  At  the  critical  point  9'  ~  94.6°,  E  (a)  changes 
qualitatively,  signals  another  phase  transition. 


2A  /3*  0  pf 


4 


A 

o 


2A  p*  0 

~P2  -2A  -p* 


-Pi  0  -p3  -2A 


(14) 


The  expressions  for  ph  p2  and  p3  are  lengthy  and  tabulated  in  appendix  B.  In  what  follows,  we  will  discuss  the 
energy  E(p,  a)  =  E0  +  Ex  separately  for  the  two  distinct  cases:  O'  <  90°  and  O'  >  90°. 

The  results  for  E  (ip)  from  the  spin  wave  analysis  are  plotted  in  the  upper  panel  of  figure  4  for  several  values 
of  O'  corresponding  to  the  Neel  ordered  phase.  One  notices  that  the  fluctuations  do  not  change  qualitatively  the 
mean  field  ground  state.  E  (p)  reaches  minima  still  at  p  =  0  or  7r  for  small  O' .  However,  as  O'  is  increased,  the 
energy  E  (p)  becomes  flatter.  Eventually,  as  O'  =  0'c~  75.0°  (the  top  curve  of  figure  4),  the  energies  for 
p  =  7r/4,  37t/4  with  proper  choice  of  a  become  degenerate  with  those  for  p  =  0,  tt.  This  signals  the 
destabilization  of  the  Neel  order  by  quantum  fluctuations.  This  occurs  around  Qfc,  before  the  Kitaev  point  is 
approached.  Note  that  in  figure  4,  only  the  region  p  G  [0,  7r/4]  (J  [37t/4,  tt]  is  shown.  Outside  this  region  (and 
also  for  o'  >  o'y  the  lowest  order  spin  wave  theory  based  on  the  Ising-like  antiferromagnetic  order  becomes  ill 
defined. 

For#7  >  90°,  the  classical  ground  state  is  continuously  degenerate  with  p  =  7r/2  but  arbitrary  ce  G  [0,  2tt]. 
Quantum  fluctuations  lift  the  degeneracy  and  select  a  long-range  ordered  ground  state  via  the  ‘order  by  disorder 
mechanism’.  Such  mechanism  for  the  special  case  of  O'  =  120°,  i.e.  the  quantum  120°  model,  has  been 
discussed  before  in  [6,  7, 34].  As  shown  in  the  lower  panel  of  figure  4,  the  same  physical  picture  continues  to  hold 
for  the  tripod  model  for  O'  ^  120°:  the  energy  E  is  minimized  at  an  =  mr/3  with  integer  n.  However,  E  (a) 
becomes  increasingly  flat  as  O'  is  decreased.  At  the  critical  point  01  =  O'i  ~  94.6°,  additional  minima  of  E  appear 
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at  a  =  an  +  tt/6.  This  indicates  that  the  long- range  order  gets  destroyed  and  replaced  by  a  new  phase 
for  0r  <  0'd. 

We  emphasize  that  the  simple  version  of  spin  wave  theory  outlined  above  is  only  intended  to  estimate  the 
lower  and  upper  critical  points  for  the  SL  phase,  0'c  and  0'd.  It  can  be  further  improved  to  properly  treat  general 
classical  spin  configurations.  We  will  not  do  it  here  because  the  large  S  expansion  by  itself  cannot  unambiguously 
determine  the  order  for  our  model  of  S  =  1/2  in  the  region  0'  >  90°. 

5.  The  dimer  phase 

Previous  theoretical  studies  of  the  quantum  120°  model  on  the  honeycomb  lattice  gave  conflicted  results.  The 
spin  wave  analysis  of  [6]  assumed  a  homogenous  ground  state  Sr  =  S0  and  found  quantum  fluctuations  prefer 
an  =  mr/3.  This  led  the  authors  to  suggest  that  the  ground  state  maybe  a  simple  ferromagnet  of  S  with  any 
choice  of  an  (i.e.  antiferromagnetic  order  in  terms  of  the  original  spin  T  or  r).  Reference  [7]  considered  more 
general  (inhomogeneous)  classical  ground  states  and  discovered  that,  within  spin  wave  theory,  the 
ferromagnetic  state  is  energetically  less  competitive  than  a c fully  packed  unoriented  loop  configuration’  with  the 
same  an  values.  The  reason  is  quite  subtle  but  argued  to  be  physically  robust:  the  loop  configuration  hosts 
maximum  number  of  zero  modes.  This  result  obtained  from  semiclassical  large-5  expansion  was  conjectured  to 
survive  in  the  limit  of  5  =  1/2,  i.e.  the  quantum  120°  model  has  a  ground  state  with  the  six-site  plaquette  order 
[7].  However,  no  evidence  of  long-range  order  was  found  in  ED  studies  where  the  spin  correlation  functions 
were  computed  for  finite  size  clusters  with  periodic  boundary  conditions  [34].  Instead,  the  ED  results  supported 
a  trial  wave  function  similar  in  spirit  to  the  short-range  resonating  valence  bond  state,  i.e.,  a  liquid  state  with 
linear  superposition  of  dimer  covering  of  the  lattice.  Therefore,  the  true  ground  state  of  the  quantum  120°  model 
was  not  settled. 

Compared  to  these  previous  works,  the  numerical  tensor  network  algorithm  used  here  takes  into  account 
quantum  fluctuations  beyond  the  lowest  order  spin  wave  theory,  works  directly  in  the  thermodynamic  limit,  and 
starts  with  unbiased  (random)  choice  of  tensors  as  variational  parameters.  It  is  capable  of  describing  both  the 
long-range  ordered  and  the  spin-disordered  ground  states.  We  find  the  ground  state  of  the  quantum  120°  model 
is  the  dimer  phase  illustrated  in  figure  5(a)  where  the  arrows  denote  the  direction  ft  of  spin  average  (r)  on  each 
site,  and  the  numbers  indicate  the  bond  energies  in  unit  of  /.  The  long-range  spin  order  we  observed  agrees  with 
the  conjecture  based  on  physical  insights  in  [7]. 

We  prefer  the  shorter,  more  descriptive  name  of ‘dimer  phase’  adopted  here  because  it  indicates  a  solid 
(crystal)  order  of ‘dimers’,  i.e.  antiferromagnetically  aligned  spins  along  the  bond  direction,  on  a  subset  of  the 
bonds.  We  propose  to  describe  the  long-range  order  using  the  vorticity  or  the  winding  number  of  the  spin 
configuration  around  each  hexagon 


v  = 


1 

•  «;+i» 

5  M 


(15) 


where;  labels  the  six  sites  of  the  hexagon.  For  example,  hexagons  marked  by  a  dot  in  the  center  in  figure  5(a) 
correspond  to  v  =  1  where  all  spins  on  the  vertices  point  radially  outwards  (corresponding  to  the  ‘loop’  in  [7]). 
The  rest  of  the  hexagons,  each  marked  by  a  cross  at  the  center,  have  v  —  — 1/2.  It  then  becomes  apparent  that 
the  hexagons  marked  by  dots  form  a  triangular  lattice  of  spin  vortices.  And  within  one  unit  cell  of  the  triangular 
lattice,  the  total  vorticity  is  zero.  Note  that  the  state  shown  in  figure  5(a)  is  energetically  degenerate  with  a  state 
where  all  the  spins  are  flipped. 

By  embedding  the  quantum  120°  model  into  the  more  general  tripod  model,  we  are  able  to  monitor  the 
suppression  of  the  dimer  order  and  its  eventual  transition  into  the  gapless  SL  phase  (phase  B)  of  the  Kitaev 
model.  The  results  are  summarized  in  figure  6.  We  observe  that  the  in- plane  magnetization  02  decreases 
continuously  to  zero  as  0'  is  reduced.  Meanwhile,  the  ground  state  energy  E  steadily  rises,  indicating  an  increased 
degree  of  frustration  as  the  Kitaev  point  is  approached.  One  can  measure  the  dimer  order  by  introducing  the 
energy  difference  6E  between  the  averages  of  two  types  of  bonds:  the  ‘happy’  bonds  (dimers)  with  antiparallel 
spins  and  the  frustrated  bonds  where  the  two  spins  form  an  angle  of  60°.  Therefore,  the  dimensionless  parameter 

rj  =  6E/E  (16) 

can  also  serve  as  the  order  parameter  for  the  dimer  phase.  As  plotted  in  figure  6,  rj  continuously  drops  to  zero  as 
6'  is  reduced  from  120°  to  95°.  Once  inside  the  SL  phase,  both  02  and  ^vanish,  and  the  bond  energies  become 
approximately  the  same  (see  figure  5(b)).  One  can  view  the  transition  from  the  SL  to  the  dimer  phase  as 
condensation  of  spin  vortices.  Equivalently,  when  0'  is  reduced,  one  can  view  the  demise  of  the  dimer  order  as 
the  melting  of  the  spin  vortex  lattice.  Note  the  bond  energy  shown  in  figure  5  features  small  fluctuations  and  does 
not  strictly  obey  C6  rotation  symmetry.  In  our  tensor  network  calculations,  no  spatially  symmetry  is  enforced  on 
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Figure  5.  Configurations  for  (a)  the  dimer  phase  (O'  =  120°)  and  (b)  the  spin  liquid  phase  (O'  =  90°)  of  the  tripod  model.  The 
number  associated  with  each  bond  is  ( Tr7  T?+e  ),  i.e.  the  bond  energy  in  unit  of /,  from  HOTRG  calculation.  The  arrows  in  (a)  depict 
the  spin  direction  on  each  site. 


the  tensors,  and  the  expectation  values  of  operators  are  computed  approximately.  The  fluctuations  are  expected 
to  decrease  as  the  bond  dimension  is  increased. 


6.  Potential  realization 

The  tripod  model  can  be  realized  from  the  following  Hubbard  model  at  half  filling  in  the  Mott  limit 


Hhub  =  -  E  C'flNe, 


^E 


(17) 


l,  7,<T<T 
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Figure  6.  The  order  parameters  of  the  dimer  phase,  02  and  rj  (defined  in  the  main  text),  and  the  ground  state  energy  £  as  functions  of 
6' .  The  vanishing  of  02  and  r\  marks  the  transition  from  the  dimer  phase  to  the  spin  liquid  phase. 


where  a  annihilates  a  fermion  with  spin  a  at  site  i.  The  direction  and  spin  dependent  hopping  matrix  is  related 
to  T7  defined  in  equation  (2)  simply  by 

£,  =  f[l/2  +  T7(0)W.  (18) 


Explicitly,  they  are  given  by 


hi  —  “(1  —  C1C2)>  £|T  =  ~0  +  clc2)> 

t  .  .  .  t  .  .  . 

^Tf  =  -fe52  -  l*i),  hi  =  ~  (^l52  +  15l), 


where  we  have  suppressed  the  superscript  7  for  brevity,  and 

Si  =  sin  6 ,  Ci  =  cos  6 ,  S2  =  sin  0  ,  c2  =  cos  (f>  . 

In  the  limit  of  U  t,  using  second- order  perturbation  theory,  we  obtain  the  following  effective  Hamiltonian 

for  Hhub 

He ff  =  jJ2Tr  (0)  V+e  (9)  -  L  (19) 

r,7  4 

which  is  nothing  but  the  tripod  model  H  ( 0 ),  up  to  a  constant  term,  with  the  superexchange  /  =  2  t2/U .  Note 
that  the  derivation  of  the  effective  compass  Hamiltonian  above  does  not  depend  on  the  details  of  the 
parameterization  of  T7  in  terms  of  6  and  <j)  .  It  follows  that  a  large  class  of  compass  models,  not  limited  to  the 
tripod  model  proposed  here,  can  be  engineered  on  honeycomb  lattice  following  the  recipe  above. 

Duan  et  al  previously  showed  that  the  Kitaev  model  can  be  realized  using  cold  atoms  on  a  hexagonal  optical 
lattice  with  extra  laser  beams  [36].  Generalization  of  their  idea  to  the  case  of  the  tripod  model  (and  other 
compass  models)  requires  spin- dependent  hopping  tGG>  controlled  by  a  non-Abelian  gauge  field  or  generalized 
spin-orbit  coupling.  Schemes  to  realize  spin-orbit  coupling  was  proposed  in  various  approaches  [37-40].  The 
realization  of  many  have  been  demonstrated  successfully  in  cold  atoms  experiments  [41—44].  For  example,  spin- 
dependent  optical  lattices  have  been  engineered  using  magnetic  gradient  modulation  [44-46].  It  seems  possible, 
but  challenging,  to  make  tGGf  spatially  dependent.  Alternatively,  the  tripod  model  proposed  here  maybe 
emulated  using  other  artificial  quantum  systems  such  as  superconducting  quantum  circuits  [47]. 


7.  Outlook 

The  tripod  model  introduced  in  this  paper  encompasses  three  well  known  models  of  quantum  magnetism:  the 
Ising  model,  the  Kitaev  model  and  the  120°  model.  We  established  its  (zero  temperature)  phase  diagram  using 
tensor  network  algorithms.  This  amounts  to  solving  a  continuum  of  frustrated  spin  models  with  spatially 
dependent  exchange  interactions.  In  particular,  we  found  an  extended  SL  phase  around  the  Kitaev  point,  and  a 
dimer  phase  for  large  values  of  angle  6r  including  the  quantum  120°  model.  The  two  quantum  critical  points 
obtained  from  tensor  network  states  agree  roughly  with  estimations  from  spin  wave  theory. 

Our  work  only  scratches  the  surface  of  the  rich  physics  contained  in  the  tripod  model.  Here  we  mention  just 
a  few  open  questions  to  be  addressed  in  future  work.  First  of  all,  it  is  desirable  to  develop  a  fieldtheoretical 
description  of  the  continuous  phase  transition  between  the  SF  phase  and  the  long-range  ordered  dimer  phase, 
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based  on  the  intuitive  picture  of  spin  vortices  introduced  in  section  5.  Secondly,  the  tripod  model,  like  other 
compass  models,  has  very  interesting  emergent  symmetry  properties  including  intermediate  symmetries 
midway  between  the  global  and  local  symmetries  [1].  Consequently,  the  excitation  spectrum  is  expected  to 
contain  zero  modes  and/  or  flat  bands.  It  is  therefore  valuable  to  understand  the  excitation  spectra  of  the  long- 
range  order  phases  by  going  beyond  the  ground  state  analysis  here.  Thirdly,  the  finite  temperature  properties  of 
the  tripod  model  deserve  a  separate  study.  The  classical  limit  of  the  tripod  model  is  known  to  be  highly 
nontrivial.  The  effects  of  thermal  fluctuations  and  the  ‘order  by  disorder’  mechanism  have  been  investigated  in 
[34]  for  the  classical  120°  model.  Finally,  we  have  only  focused  on  the  case  of  /7  =  /  here.  From  the  Kitaev 
model,  we  know  that  a  gapped  SL  phase  (phase  A)  takes  over  when  the  asymmetry  in  J1  grows  large.  Thus  one 
expects  that  further  generalization  of  the  tripod  model  to  general  values  of  J1  may  uncover  new  interesting 
phases. 

To  conclude,  we  hope  our  results  can  stimulate  further  application  of  tensor  network  algorithms  to 
frustrated  spin  models  as  well  as  spin-orbital  models  describing  transition  metal  oxides.  We  also  hope  our 
introduction  of  the  tripod  model  can  inspire  alternative  proposals  to  extend  the  Kitaev  model  or  realize  compass 
models  in  artificial  quantum  systems  such  as  cold  atoms  on  optical  lattices  or  superconducting  circuits. 
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Appendix  A.  Tensor  network  algorithms 

The  finite  PEPS  algorithm  is  a  powerful  numerical  approach  for  two-dimensional  quantum  spin  systems  [25- 
27].  Lor  the  tripod  model,  we  construct  the  usual  PEPS  wave  function  starting  from  six  random  rank- four 
tensors  with  virtual  bond  dimension  D  =  3.  The  tensors  are  then  optimized  through  recursive  imaginary- time 
evolution  with  time  step  r  =  0.01  on  all  the  links.  Once  the  wave  function  is  converged,  we  calculate  the  ground 
state  energy  and  the  expectation  value  of  the  combined  order  parameter  O  =  yjo?  +  022.  The  results  are 
shown  in  figure  7  for  three  typical  values  of  0'  (corresponding  to  the  three  different  phases  found  in  the 
thermodynamic  limit).  Lor  the  small  system  size  considered  here  (six  sites  with  periodic  boundary  conditions), 

O  vanishes  as  the  wave  function  converges  to  the  ground  state.  Nonetheless,  a  noticeable  peak  of  O  during  the 
time  evolution  can  serve  as  the  indication  for  spin  order.  Lor  all  the  three  cases,  the  energies  converge  to  the  ED 
result  up  to  a  relative  error  of  10-3  (the  upper  panel  of  figure  7).  The  peaks  of  O  in  the  lower  panel  of  figure  7  for 
the  cases  0'  =  70°  and  0'  =  120°  point  to  the  Neel  order  phase  and  the  dimer  phase  respectively,  while  the 
monotonic  decay  of  O  for  the  case  0'  =  90°  suggests  a  SL  state. 

Within  the  many  variants  of  tensor  network  algorithms,  a  typical  way  to  find  the  phase  diagram  of  quantum 
spin  models  is  the  infinite  PEPS  (iPEPS)  method  [28, 29].  The  iPEPS  ansatz  on  the  honeycomb  lattice  usually 
proceeds  by  mapping  the  lattice  to  a  square  lattice  and  evaluating  the  effective  environment  by  contraction 
schemes  such  as  infinite  matrix  product  states  [28]  or  corner  transfer  matrices  [48, 49].  For  instance,  the  phases 
of  Kitaev-Heisenberg  model  [50]  and  the  SU(4)  symmetric  Kugel-Khomskii  model  [51]  have  been  studied  via 
the  iPEPS  ansatz  with  a2  x  2  or  4  x  4  unit  cell.  However,  the  contraction  scheme  for  a  six- site  (hexagonal)  unit 
cell  on  the  honeycomb  lattice  is  tedious  and  expensive,  especially  for  the  corner  transfer  matrices  scheme. 

For  this  reason,  we  adopt  the  simple  tensor  update  scheme  and  evaluate  the  contraction  using  the  HOTRG 
method  as  explained  in  the  main  text.  The  simple  update  [32]  generalizes  the  TEBD  [33]  technique  to  two- 
dimensional  quantum  systems  by  introducing  the  bond  vectors  to  represent  the  mean- field  environment  for 
local  tensors.  We  set  the  imaginary  time  step  r  =  0.01  and  the  number  of  iterations  is  generally  around  105 
(smaller  time  step  does  not  improve  the  numerical  result  significantly).  The  accuracy  of  the  HOTRG  method  is 
controlled  by  the  virtual  bond  dimension  D.  By  systematically  increasing  D,  the  quantum  entanglement  between 
neighboring  sites  is  better  taken  into  account,  yielding  a  more  accurate  ground  state.  For  example,  figure  8  shows 
that  the  order  parameter  02  vanishes  when  D  is  increased  to  8  in  the  region  0r  <  94°,  suggesting  a  SL  ground 
state.  One  notices  that  the  variations  of  the  ground  state  energy  with  O'  within  the  SL  phase  is  larger  than  those  in 
the  long-range  ordered  phase,  especially  for  smaller  D  values.  This  is  due  to  the  strong  quantum  fluctuations 


11 


IOP  Publishing 


New  J.  Phys.  18  (2016)  053040 


HZou  etal 


N 

Figure  7.  Finite  PEPS  results  for  a  six-site  cluster  showing  the  errors  of  ground  state  energy  (relative  to  that  from  the  exact 
diagonalization)  and  the  order  parameter  O  for  three  different  values  of  O' . 


Figure  8.  The  ground  state  energy  E  and  the  order  parameter  02  computed  from  HOTRG  at  virtual  bond  dimension  D  =  4,  6,  8and 
truncation  dimension  x  =  8. 


intrinsic  to  the  SL.  Cross-checking  the  HOTRG  calculations  here  to  those  using  Second  Renormalization  Group 
[52]  which  takes  into  account  the  entanglement  between  the  system  and  the  environment  deserves  a  future 
study. 
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Appendix  B.  Spin  wave  theory 


For  the  two  long-range  ordered  phases,  the  1/5  corrections  to  the  energy  are  obtained  by  diagonalizing  the 
matrix  equation  (14)  with  different  form  factors  /3h  (32,  and  /?3.  For  the  Neel  ordered  phase,  they  are  given  by 

3 

A  =  EA;  +  ifl;)2eik^, 

j= 1 

3 

A  =  E  A;  +  ia;)2e_ikA 

i= 1 

3 

A  =  -  E Aj  +  ^2)eiH- 

;'=i 


And  for  the  dimer  phase 

3 

A  =  E  A;  +  i^;)2eik'eh 

i=i 

3 

A  =  E  A;  +  ibj)2e~ik^, 

j= 1 

A  =  -  E^2  +  bj)eik'^. 

j=l 


Here,  a;-  and  bj  are  related  to  the  parameter  0 ,  (p,  and  ce  defined  in  the  main  text  through 


a\  =  cos  6  cos  (p  cos  ce  —  sin  6  sin  (p, 


a2  =  cos  0  cos  ip  cos 

a3  =  cos  6  cos  <p  cos 
bi  =  —  cos  6  sin  a, 
b2  =  —  cos  6  sin 


(“-tH 

DtH 


sin  #  sin  <p, 


sin  0  sin  <p, 


cos  6  sin 


in(a  -  ?)• 
in(“  -  f  )■ 
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